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Liquids spreading over a solid substrate under the action of various forces are known to exhibit a 
long wavelength contact line instability. We use an example of thermally driven spreading on a 
horizontal surface to study how the stability of the flow can be altered, or patterns selected, using 
feedback control. We show that thermal perturbations of certain spatial structure imposed behind 
the contact line and proportional to the deviation of the contact line from its mean position can 
completely suppress the instability. Due to the presence of mean flow and a spatially nonuniform 
nature of spreading liquid films the dynamics of disturbances is governed by a nonnormal evolution 
operator, opening up a possibility of transient amplification and nonlinear instabilities. We show that 
in the case of thermal driving the nonnormality can be significant, especially for small wavenumber 
disturbances, and trace the origin of transient amplification to a close alignment of a large group of 
eigenfunctions of the evolution operator. However, for values of noise likely to occur in experiments 
we find that the transient amplification is not sufficiently strong to either change the predictions of 
the linear stability analysis or invalidate the proposed control approach. 
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I. INTRODUCTION 

Driven spreading of liquid films is a process which oc- 
curs in numerous industrial coating applications, such 
as spin-coating of hard drives, so understanding its dy- 
namics and learning to control it is very important. For 
instance, instabilities which arise during the spreading 
of the liquid on the solid substrate can lead to nonuni- 
form coverage, adversely affecting the quality of produced 
coating. Driven spreading of thin films and patterning 
also have important implications for microfluidics. 

Driven spreading of liquid films under the action 
of gravity^'^, centrifugal acceleration"^, thermocapillary 
effects'', or combination thereof^ has been extensively 
studied in the literature. Stability analysis of such flows 
has attracted the most attention and the mechanism of 
the linear instability is now well understood. Consider- 
able progress has also been reached in feedback control 
of flat liquid layers^"*, whose dynamics is governed by 
normal differential operators. The attempts to influence 
the stability of spreading films have so far been limited 
to non-feedback control achieved through either impos- 
ing an externally generated counterflow^ or chemically 
patterning the substrate^°'^^. 

This study represents the first theoretical treatment of 
the feedback control problem for spreading films. The 
spatially nonuniform nature of spreading films and the 
presence of mean flow make the control problem much 
more difficult compared to the case of flat stationary 
films, because the dynamics of disturbances in spread- 
ing films is governed by a nonnormal evolution operator 
and thus requires a more careful analysis. For instance, 
liquid films flowing down an incline have been found to 
develop a contact line instability when the linear analysis 
predicted stable evolution^^. We derive the slip model of 



thermally driven spreading and use it to show that the 
contact line instability can be suppressed using adaptive 
thermal perturbations which depend on the distortion of 
the contact line. (This type of feedback is chosen because 
it is easiest to implement experimentally with sufficient 
spatial and temporal resolution via optical means^'^.) Al- 
though the results of the following analysis should be ap- 
plicable regardless of the driving force, we concentrate 
our attention on the case of thermal driving. 

The layout of the paper is as follows. The slip model of 
thermally driven liquid films is derived in Section II and 
its linear stability analysis is conducted in Section III. 
The validity of linear stability analysis and transient ef- 
fects are considered in Section IV. Section V presents the 
proposed algorithm for feedback control of the contact 
line instability and Section VI contains the conclusions. 



II. SLIP MODEL FOR THERMAL SPREADING 

We consider the spreading of a thin layer of partially 
wetting liquid on a horizontal substrate (see Fig. 1). The 
spreading process is conventionally described using the 
lubrication approximation^^, with the horizontal velocity 
governed by the Stokes equation 



ti-dzz^ = Vp, 



(1) 



where fi is the dynamic viscosity, p is the modified pres- 
sure, V = {dx, dy) is the 2-dimensional gradient operator, 
and the vertical velocity is neglected. 

It is well known^^'^^ that the standard no-slip bound- 
ary condition at the liquid-solid interface results in a 
stress singularity at the contact line. The only approach 
explored in the literature for the thermally driven case 
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FIG. 1. Spreading liquid film on a solid substrate. 



was to relieve this singularity by introducing a thin pre- 
cursor film^. However, as we intend to use the position 
of the contact line in our control algorithm later, the pre- 
cursor model becomes inconvenient. Instead we choose 
to specify a microscopic contact angle and relieve the 
stress singularity by employing a partial slip boundary 
condition 



3h 



(2) 



at the bottom of the liquid layer, where h is the local 
thickness of the film, and a is the phcnomcnological slip 
coeflScient. The slip boundary condition (2) was orig- 
inally introduced by Greenspan^^ for modeling the un- 
forced spreading of liquid drops and later used to model 
the forced spreading of liquid films under the action of 
gravity^^. At the free surface the standard stress balance 
boundary condition 



(3) 



applies, where a is the surface tension coefficient. Solving 
(1) subject to these boundary conditions we obtain the 
horizontal velocity 
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In order to make the phenomenological boundary condi- 
tion (2) consistent with the physics of the flow, in (4) 
we have dropped an unphysical term aVa/S/i/i which 
diverges for vanishing film thickness. This divergence is 
not necessarily a cause for alarm as the continuous ap- 
proximation underlying the Stokes equation itself breaks 
down in this limit. The shear stress also becomes poorly 
defined for very thin films. Our choice of the functional 
form for the horizontal velocity (4), therefore, amounts 
to picking an appropriate phenomenological model in the 
region where the continuous description of the fiow be- 
comes invalid and fluctuations become important. 

For a film which is sufficiently thin the hydrostatic 
pressure can be ignored, so that the modified pressure 
is given by the normal component of the surface tension 
p = —an = —aV^h. Substituting (4) into the mass con- 
servation condition 
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(V • v)dz 



(5) 



and integrating we obtain an evolution equation for the 
thickness: 



dth = -V • 



(6) 



Now consider the situation which arises when the sub- 
strate covered by the liquid film is subjected to a linear 
temperature gradient in the x-direction. Assuming that 
the surface tension changes linearly with temperature 6, 



cr{x) = (7{0a) + xdxOdgcr = cto 



(7) 



and neglecting the variation in a in the second term of 
(6), which produces subdominant contribution (see, e.g., 
the discussion in Ref. 5), we obtain 

dth = ^9xh^ - ^dx[{ah -\- h^){dxxxh + d^yyh)] 

~ + h^){dxxyh + dyyyh)]. (8) 

Ojjj 

We can absorb most parameters into the spatial and 
temporal scales by introducing the nondimensional vari- 
ables t' = t/T, x' = x/X, y' = y/X, and h' = h/H. 
The vertical length scale H is defined by the charac- 
teristic thickness of the film and sets both the horizon- 
tal length scale X = [2cro-ff^/3r]^/^ and the time scale 
T = 2^X/tH. After defining the dimensionless slip co- 
efficient a' = a/H^ and dropping the primes we obtain 



dth = dxh^ - dx[{ah + h^){dxxxh + dxyyh)] 

- dy[{ah + h^)(dxxyh + dyyyh)]. 



(9) 



The obtained equation has the same form as the one de- 
scribing gravity driven rather than temperature driven 
films (see, e.g., equation (33) in Ref. 17), with the excep- 
tion that h^ in the first term on the right-hand-side is 
replaced with ah + h^. This similarity of the structure of 
the governing equations suggests that the feedback con- 
trol problem can be treated in the same way regardless 
of the nature of the driving force. 

The liquid spreads in the direction opposite to the tem- 
perature gradient, so the motion of the contact line is 
most conveniently described in the reference frame mov- 
ing with speed u towards negative x. In this frame 
the evolution equation possesses a transversely uniform 
steady state solution, which gives the asymptotic film 
profile for constant flux boundary conditions. Substitut- 
ing h{x,y,t) = ho{x + ut) into (9) and integrating once 
we obtain 



uho -hl + {aho + hDh^ = d, 



(10) 



where prime indicates the differentiation with respect to 
the X coordinate. The constants u and d can be deter- 
mined from the appropriate boundary conditions. Fol- 
lowing Spaid and Homsy-'^^ we require that the film thick- 
ness vanish at the contact line, ho{0) = 0, and spec- 
iiy the slope /io(0) = c = (X/iJ)tan0, where <p is the 
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FIG. 2. Asymptotic film profiles produced by the precursor piQ 3 Asymptotic film profiles with a fixed height of the 

model (12) with b = 0.0119 and the slip model with a = 0.1 capillary ridge hmax = 2.75 and different values of the slip 

and c = 1.53. The parameters were chosen to produce a coefficient, 
capillary ridge of height hmax = 2. 



microscopic dynamic contact angle in the dimensional 

variables. Furthermore, the constant flux boundary con- 
dition at the tail of the film far away from the contact 
line gives /lo(oo) = 1 (this choice corresponds to using 
the dimensional tail thickness as the vertical Icngthscalc, 
H) and /iq"((X)) = 0. These boundary conditions require 
d = and u= 1, and consequently 



ho-1 
hl+a 



(11) 



The sohition to this equation describes the height profile 
of the spreading film once the distance from the contact 
line to the reservoir becomes sufficiently large. 

At this point it is appropriate to make two comments 
regarding the structure of the lubrication equations pro- 
duced by the slip model. First of all, from a mathematical 
perspective, introduction of partial slip at the liquid-solid 
interface is equivalent to a regularization procedure for 
a singular differential equation. Dropping the diverging 
term in (4) is equivalent to keeping the regularization 
parameter a only in the terms in (6) responsible for the 
singular behavior of the solution at the contact line. We 
could have kept all the terms in (4) and (6) just as well. 
Neither (9) nor (11) would have changed. 

Second, (11) is very similar to the equation produced 
by the precursor model of Kataoka and Troian^ 



,„ _ {hp - l){ho - b) 



(12) 



where b is the thickness of the precursor layer. In fact, 
the solutions of the two equations are virtually indistin- 
guishable (see Fig. 2) for the proper choice of parameters 
(precursor thickness, slip coefficient, and contact angle). 
This is hardly a surprise, as the two equations become 
identical in the limit 6 — > and a — > 0, and proves that 
the final result is essentially independent of the regular- 
ization procedure used to get rid of the stress singularity 
at the contact line. 



The slip model has two independent parameters (c and 
a), while the precursor model has only one (6). A closer 
inspection shows that the gross features of the slip model 
are, in fact, determined by a single parameter (say, the 
height of the capillary ridge). As Fig. 3 shows, one can 
change the value of the slip coefficient by two orders of 
magnitude without causing noticeable changes in the film 
profile (the contact angle will also change in response to 
the changes in a to keep the height of the capillary ridge 
constant). The second parameter fine tunes the shape 
of the film in the immediate vicinity of the contact line. 
The slip model is, therefore, not only more convenient 
for the purpose of control (it explicitly defines the easy 
to measure position of the contact line), but it is to some 
extent more flexible than the precursor model. 

Finally, let us explore the dependence of the film pro- 
file on the parameters. As we have just seen, the profile 
is rather insensitive to the changes in the slip coefficient, 
so let us fix it and vary the contact angle (see Fig. 4). 
The first important observation is that regardless of the 
magnitude of the contact angle, the capillary ridge never 
disappears (this turns out to be the case for any reason- 
able value of the slip coefficient, i.e., a < 1). This is 
in contrast with the results obtained for gravity driven 
films, where the capillary ridge completely disappears at 
small inclination angles^^. Second, the height of the ridge 
is a monotonically increasing function of c and becomes 
essentially independent of the contact angle for c < 0.3. 
The dependence on the slip coefficient is more subtle: the 
value of a controls the minimum height of the capillary 
ridge which grows with decreasing slip coeSicient. 



III. CONTACT LINE INSTABILITY 

Linear stability of the asymptotic solution ho{x + ut) 
can be determined in a standard way^^'"''^^. Since this 
solution is uniform in the transverse direction, the lin- 
earized equation describing the dynamics of small dis- 
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FIG. 4. Asymptotic film profiles for different values of the 
contact angle. The slip coefficient is held fixed at a = 0.01. 



turbances can be partially diagonalized by Fourier trans- 
forming it in the ?/-direction. By substituting 

h{x, y, t) =ho{x + ut) + eg{x + ut, t)e''^'" (13) 

into (9) and retaining terms of order e, we obtain 

dtg = L{q)g, (14) 

where L{q) = Lq + q^Li + g^L2 is a 4*^ order differential 
operator defined via 

Log = -[{l- 2ho + {a + 3hl)h'^'}g + {aho + h^g'"]', 

Lig = [{a + 3hl)h'og]' + 2{aho + hl)g" , 

L2g = -{aho + hl)g. (15) 

The boundary condition on g at the contact line 

K{0)giO) = KiO)g'{0) (16) 

can be obtained by requiring the perturbed solution to 
have the same contact angle as the unperturbed solution. 
The other three boundary conditions can all be imposed 
at the tail, say, 



5(00) = 5' (00) = g"{oo) = 0. 



(17) 



Even though we cannot find the eigenvalues of L(q) for 
arbitrary q analytically, for long wavelength disturbances 
we can use perturbation theory to get the leading order 
(in q"^) terms. This requires finding the eigenfunctions of 
Lq and its adjoint, Lq. Taking the second derivative of 
(10) we obtain 



Lo h'o = 0, 



(18) 



so that go = h'^ is an eigenfunction of Lq with eigenvalue 
/3o = 0. The adjoint operator is defined via 

Llf =[l-2ho + {a + ihl)K']f' 

- [{aho + hDf]'", (19) 



so its respective eigenfunction is just a constant, say, 
/o = 1. In fact, these are generic results with deep phys- 
ical meaning. Identical relations between the asymptotic 
state and the leading eigenfunctions were obtained, e.g., 
for gravity driven films using the precursor model^^'^^. 
The relation for go is due to the fact that equations for 
the asymptotic state are translationally invariant in the 
direction of the flow (this reflects an arbitrary choice in 
the position of the contact line) , while the relation for /o 
is the consequence of the divergence form of (5), which 
reflects mass conservation. 

According to the perturbation theory the leading 
eigenvalue has the following g-dcpcndence: 



/3o(9) 



ir fo9odx 



+ 0(g4 



Using (11) this can be reduced to 



/3o(g) 
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ho{ho-l)dx + 0{q^). 



(20) 



(21) 



As this is the largest eigenvalue, its sign determines the 
stability of the asymptotic state. It is easy to see that, 
if the asymptotic profile were monotonic, < /lo < 1, 
the integral would be strictly negative and the system 
would be linearly stable with respect to long wavelength 
disturbances. However, the presence of a capillary ridge 
makes the integral positive, showing that the increased 
mobility of the ridge provides the mechanism for the long 
wavelength instability in the thermally driven case. This 
mechanism has been originally conjectured by Kataoka 
and Troian based on the energy analysis of the precursor 
model^. Equation (21) gives an explicit condition on the 
shape of the capillary ridge which determines the onset 
of instability, and echoes a similar result obtained for the 
case of gravity-driven flows^'^^. 

Substituting g{x,t) = /iq(x) exp[/3o(fj')i] into (13) we 
notice that for small disturbances the right hand side 
represents the first two terms of the Taylor expansion of 
a distorted asymptotic state ho{x + ^ + ut), where 



ee 



(22) 



is the deviation of the contact line from the mean. In 
fact, the marginal translational mode go is not the only 
eigenfunction of Lo- There is an infinite discrete spec- 
trum of eigenvalues /?„ and eigenfunctions gn- Therefore, 
in the presence of an arbitrary disturbance (22) should 
be replaced with 

1 

^{y,t) = -J29n{0) e„(g)e^«^+'5"(«)*d(z. (23) 

As the asymptotic state is imstable, the amplitude of 
the distortion will grow exponentially in time and even- 
tually the contact line will form equally spaced "fingers" . 
In order to calculate the wavenumber of the pattern we 
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numerically compute the eigenfunctions and eigenvalues 
of the evolution operator L{q). This is most easily ac- 
complished by discretizing the eigenvalue equation on a 
spatially nonuniform mesh to properly resolve the rapid 
change in the solutions near the contact line. As both 
the asymptotic state ho and the eigenfunctions gn expo- 
nentially flatten for large x, a truncated domain can be 
used, so that the boundary conditions (17) arc imposed 
at finite distance away from the contact line (we used 
Ix = 80 in most of the calculations). In order to com- 
pute the asymptotic state we used a shooting method on 
a nonuniform mesh with roughly 10** points. As h = 1 
is an unstable fixed point of (11), numerical integration 
(we used the 4'^ order Runge-Kutta method) cannot be 
performed beyond Ix « 30 using double precision arith- 
metics. To determine the solution for a longer interval 
we used a mcthod^^ in which the numerical solution is 
extended using the analytical solution of (11) for ft, « 1: 

h{x) = 1 + [acos{V3Kx) + bsm{V3Kx)]e~'^'' H , (24) 

where n — {1 + a)~^^^/2 and a, b are constants deter- 
mined by a least squares fit. The eigenvalues and eigen- 
functions where then computed using a built-in function 
in Matlab on a 1025-point mesh (finer resolutions did not 
change the eigenvalues by more than about 5%). 

The results of our calculations for a typical choice of 
parameters arc presented in Fig. 5. The fastest grow- 
ing disturbance is found to have a transverse wavenum- 
ber Qmax lying between 0.29 and 0.34. This wavenum- 
ber decreases with the thickness of the capillary ridge 
hmax and gives the characteristic wavelength of the fin- 
gering pattern Xmax = '^^^/Qmax, which ranges between 
18.5 and 21.7, in excellent agreement with the predictions 
of both the precursor model^ and the experiments^'^*^ 
which found the dimensionless wavelength to be between 
18 and 22. The growth rate of the most unstable mode 
increases with hmax, varying between 0.12 and 0.17 for 
the range of contact angles considered here. Brzoska et 
aL^" have obtained an experimental value of about 0.15, 
which is also consistent with the theory. As the thickness 
of the ridge has not been determined in experiments, it is 
impossible to make a more direct comparison, but never- 
theless these results can be used to establish the ranges 
of parameters relevant for experimental conditions. 

IV. LINEAR STABILITY AND TRANSIENT 
DYNAMICS 

We have determined earlier that the capillary ridge 
is present for any reasonable choice of parameters, so 
the asymptotic state of a thermally driven film is always 
linearly unstable and the contact line instability will in- 
evitably set in as the asymptotic state is approached. 
However, as a quick comparison of (15) and (19) shows, 
the evolution operator L{q) is nonnormal, and we have to 




FIG. 5. The twelve leading eigenvalues of L{q) for a = 0.01 

and (a) c = 1 (fe„a, = 1.93), (b) c = 1.95 (ft„a. = 2.75). 
In both cases /3o(0) and /3i(0) are real, while the other ten 
eigenvalues come in complex conjugate pairs. 



consider the possibility that the transient dynamics asso- 
ciated with nonnormality could change the predictions of 
the linear stability theory concerning the growth rates of 
disturbances. In fact, transient effects could be quite sig- 
nificant. For instance, turbulence in channel fiows arises 
at values of the Reynolds number well below the critical 
one predicted by the linear stability analysis^ ^. Both in 
the driven films and in channel flows nonnormality arises 
due to a significant mean flow, so it is natural to expect 
that transient behavior could be important for driven liq- 
uid films as well. Indeed, a disagreement between theo- 
retical and experimental predictions of the most unstable 
wavelength in gravity driven films at low angles of incli- 
nation has been attributed to transient dynamics^^. 

Linear stability analysis presented in the previous sec- 
tion is based on the assumption that the nonlinear terms 
are negligible at all times. If the disturbances were small 
and their dynamics were governed by a normal evolution 
operator, this assumption would have been well justi- 
fied. For instance, when the system is stable, the eigen- 
values predict both the short term and the long term 
dynamics. The situation can change dramatically when 
the evolution operator becomes nonnormal, as the eigen- 
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FIG. 6. Transient amplification 7(9, t) as a function of time 
for a = 0.01, c = 1 and different values of the wavenumber. 



values become poor predictors of the short term dy- 
namics. Inclusion of the nonlinear terms in (14) has 
two major consequences. First of all, nonlinear terms 
couple the dynamics of modes with different transverse 
wavenumbers^*. Second, nonlinear terms produce devia- 
tions from the asymptotic state which can be transiently 
amplified due to the nonnormality of the linearized evolu- 
tion operator^^. A combination of these two effects can 
lead to a nonlinear instability which can compete with 
the linear instability. 

There are two scenarios which could invalidate the re- 
sults of linear stability analysis. In the simplest (purely 
linear) scenario^^'^^, an initial disturbance with trans- 
verse wavenumber qq and magnitude ||gf(a;,0)|| = ^ could 
be transiently amplified by a factor 7(^0) to produce a 
disturbance with magnitude 7^ = 0(1). If this amplifica- 
tion occurs on a time scale shorter than 1/ Poilmax) and 
^ > 7~^, the transient effects will dominate and a distor- 
tion of the contact line with the wavenumber go rather 
than Qmax will result. 

In a more complicated scenario^ an initial distur- 
bance is transiently amplified by the linear part of the 
evolution operator, while the nonlinear terms produce 
secondary disturbances which are further transiently am- 
plified. This could lead to a positive feedback loop boot- 
strapping a nonlinear instability, provided the secondary 
disturbances contain the wavenumber of the initial dis- 
turbance and have the magnitude which is at least as 
large. It is easy to check that the nonlinear evolution 
operator will only contain terms quadratic, cubic, and 
quartic in g{x, t) in addition to the linear terms kept 
in (14). Only the cubic terms will contain the original 
wavenumber, so the secondary disturbances produced by 
the quadratic and quartic terms will not be further tran- 
siently amplified. (Initial disturbances with go = repre- 
sent the only exception, but they do not lead to distortion 
of the contact line.) An initial disturbance of magnitude 
^, which is transiently amplified by a factor 7, will pro- 
duce a secondary disturbance of magnitude 0((7C)'') via 
the cubic nonlinearities. The secondary disturbance will 
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FIG. 7. Dependence of the ma^Kimal transient amplification 
7(g) on the wavenumber for a = 0.01 and different values of 
the contact angle. 



exceed the primary disturbance only when (7^)'^ ^ ^) so 
that a self-sustaining nonlinear instability becomes pos- 
sible for ^ > 7~^/^. 

The maximal transient amplification in the stable band 
can be defined in a conventional manner^^: 



, Ux,t)\\ 
7(g, t) = sup 

g{xfi) \\9{xM\ 



-Lt] 



(25) 



In the unstable band the fast transient amplification will 
be followed by the slower exponential growth. Factor- 
ing out the exponential growth one obtains the following 
upper bound: 



7(9, i) 



sup 

g{x,0) 



-Pot 



(26) 



Numerical calculations show that in the unstable band 
j{q,t) is a monotonically increasing function of time, so 
the maximum is reached for t 00. The time depen- 
dence for several values of the transverse wavenumber is 
shown in Fig. 6. As L — 0o and L have the same eigen- 
functions, we can easily calculate the matrix elements of 
the operator U (t) = exp[(X — f3o)t]: 



Uynn (^) 



Jo 



(27) 



As only the clement with m ~ n = survives for large 
times, the maximal transient amplification is achieved 
for the "optimal" initial disturbances equal to multiples 
of fo- The evolution amplifies these disturbances and 
transforms them into multiples of the leading eigenfunc- 
tion go. Exponential growth with rates predicted by the 
linear stability analysis sets in rather quickly as the time 
scales of the exponential and transient growth arc of the 
same order of magnitude (roughly between 2 and 9 nondi- 
mensional units, depending on the wavenumber). 

The ultimate test of the importance of nonlinear terms 
and transient dynamics is provided by a direct calculation 
of the transient amplification factor 
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FIG. 8. The eigenfunctions go{x) through g2i{x) for q = 0, 
a = 0.01 and c = 1. Only a portion of these eigenfunctions, 
computed on the interval of length Ix = 80, is shown. The 
eigenfunctions are normalized such that \gn\'^dx = 1. 

7((7) = max7(g,t). (28) 

The numerical results are presented in Fig. 7. For a fixed 
contact angle, 7(g) is the largest for zero wavenumber dis- 
turbances (which do not lead to distortion of the contact 
line) and quickly decreases with increasing q. The max- 
imum transient amplification increases with the contact 
angle and can become quite significant for typical exper- 
imental parameters (7(0) ~ 70 for c = 1, 7(0) « 100 for 
c = 1.95). However, the effective value for a finite sys- 
tem will likely be in the range of a few tens. The strong 
transient amplification at q = can be easily traced to a 
very close alignment of a large group of eigenfunctions (72 
through g2i. As Fig. 8 shows, their shapes are extremely 
similar. It is, therefore, appropriate to associate the tran- 
sient behavior with a whole group of stable eigenfunc- 
tions. The size of the group increases with c, increasing 
the degree of nonnormality. In contrast, for channel flows 
apparently only a couple of near-marginal eigenfunctions 
become closely aligncd^^. 

The transient amplification at small q was found to de- 
pend rather sensitively on the length of the domain used 
in the numerical calculations. For instance, the values of 
7 obtained for — 30 were generally about a half of those 
obtained for = 80. At about Ix = 80 the dependence 
levelled off, and further increase in resulted in large 
fluctuations in the eigenvalues due to numerical inaccu- 
racies resulting from strong nonnormality of the matrix 
produced by discretizing L{q). This sensitivity can be ex- 
plained by the shape of the adjoint eigenfunctions. While 
the eigenfunctions gn{x) decay exponentially for large x, 
their adjoints fnix) ^row exponentially (see Fig. 9). As a 
result, small inaccuracies in the boundary conditions at 
the tail of the film introduced by truncating the computa- 
tional domain can have a significant effect on the adjoint 
eigenfunctions, thus affecting the transient amplification, 
which depends on both gn{x) and fn{x). 

The minimal noise level required to trigger the nonlin- 



FIG. 9. The adjoint eigenfunctions fo{x) through fii{x) for 
q — 0, a = 0.01 and c — 1. The loading oigenfunction fo{x) is 
constant, aside from a small region near x = Ix, where it ad- 
justs to the Dirichlet boundary condition. The eigenfunctions 
are normalized such that f°° fngmdx = Snm- 

ear instability according to the second scenario is about 
0.1% of the total film thickness for 7 = 100 (the first sce- 
nario requires an even higher level of noise). A more ac- 
curate calculation of the noise threshold is likely to raise 
this level much higher (to maybe a few percent) by taking 
into account the fact that (26) determines the maximal 
amplification achieved for a specially chosen initial con- 
dition, while the secondary disturbances will generically 
be amplified less strongly. In typical experimental condi- 
tions the noise is likely to be substantially smaller than 
one percent of the film thickness. The existing experi- 
mental data^'^° agrees with the predictions of the linear 
theory rather well, supporting the conclusion that the 
transient effects are relatively weak and, therefore, the 
modal linear stability analysis of the previous section ac- 
curately describes the dynamics. 

V. FEEDBACK CONTROL OF THE CONTACT 
LINE INSTABILITY 

Now that we understand the limits of the linear stabil- 
ity analysis let us consider the control problem. Can the 
contact line instability be suppressed, or alternatively, 
can a pattern with a desired wavelength be imposed 
by applying feedback? In principle, the answer seems 
to be clear, as feasibility of feedback control of several 
other types of instability (buoyancy^, thcrmocapillary^, 
evaporative®) in liquid layers has been demonstrated. Al- 
though one might hope that the developed control meth- 
ods could be adapted for suppressing the contact line 
instability, in reality the spreading films turn out to be 
dramatically different. 

The existing control methods have been developed for 
stabilizing flat films with no mean flow, i.e., steady, uni- 
form target states. The evolution operators describing 
the dynamics of small disturbances about such states 
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have eigcnfunctions whose horizontal dependence is given 
by Fourier modes, which are normal to each other for any 
realistic boundary conditions (in principle, one can design 
the boundary conditions such that the eigenfunctions will 
not be normal even in this case^'*). Not only does it mean 
that all horizontal modes become uncoupled, so any mode 
can be controlled independently of the others, there is no 
nonnormality, so there are no transients and the linear 
stability analysis is unconditionally valid. 

The target state ho{x + ut) in the present problem is 
nonuniform in the direction of the flow. As a result, the 
differential operator L does not fully diagonalize and the 
control problem becomes vastly more complicated. Not 
only are all the modes in the system coupled, feedback 
applied to suppress one mode generally affects all other 
modes, so an infinite-dimensional problem has to be con- 
sidered from the outset. However, even if these problems 
are resolved and a feedback making the dynamics asymp- 
totically stable is found, there is no guarantee that the 
transient effects will not invalidate the whole analysis. 

Let us repeat the linear stability analysis of the spread- 
ing film, but now in the presence of feedback. First we 
make use of the simplification afforded by the uniformity 
of the target state in the transverse direction, which al- 
lows us to partially diagonalize the evolution operator. 
For the moment we restrict our attention to monochro- 
matic disturbances eg{x + ut, t) exp{iqy). Since the fiow 
is driven by the gradient in the temperature (and hence 
surface tension), the stability of the fiow is most easily 
altered by varying the temperature field behind the con- 
tact line. Suppose we modify the temperature profile by 
adding a perturbation 

Ae{x, y, t) = -€T{dea)-h{t)w{x + ut)e"^y, (29) 

where the transverse wavelength q is the same as that 
of the disturbance and s{t) and 'w{x) are some functions 
determining the temporal and spatial profile of the per- 
turbation, which will be determined later. Consequently, 
(7) and hence (14) will be modified to account for the 
variation in the surface tension transversely to, as well 
as along, the direction of the fiow. At order e instead of 
(14) we obtain 

dtg = Log + s[Ni{hn)w + N2{ho)w']' 

+ q^[Lig - sN2{ho)w] + q^L^g, (30) 

where the influence functions 

Ni{ho) = \{hl-ho), 

N2{ho) = hl+'^{aho + hl)h'i (31) 

determine the effect of the imposed thermal perturbation. 
To get a sense of the dynamics of different modes in 

the presence of feedback, wc expand the disturbance in 
the basis formed by the eigenfunctions of Lq, 



g{x,t) = Y.^rn{t)gm{x), (32) 

m 

and make the strength of the applied perturbation pro- 
portional to the magnitude of the distortion of the con- 
tact line (with a proportionality constant fc, called the 
gain, to be determined later), 

sit) = fc^M ^ Grn{t)gm{Q). (33) 
c c ^-^ 

m 

Multiplying (30) by /* and integrating from to oo we 
obtain an infinite system of ODEs describing the dynam- 
ics of individual modes: 

= P^G„ + '^^{Anm + Bnm + q^Cnm)Gm-, (34) 

m 

where (assuming that all adjoint eigenfunctions are nor- 
malized such that /(,°° fnOrndx = 5nm) 

Anm = r f:[m{ho)w + N2{ho)wydx, 

c Jo 

B„m= f f*Ligmdx - fc ^""^^^ f f*N2{ho)wdx, 
Jo c Jo 

/•oo 

Cnm = / fnL2gmdx. (35) 

Jo 

As Fig. 5 shows, the uncontrolled system possesses a 
single unstable eigenvalue well separated from the rest of 
the spectrum. One can, therefore, expect that the dy- 
namics of small disturbances should be well described by 
a single mode truncation of (34) in the absence of feed- 
back. The same is not generally true when the feedback 
is applied. As all modes are coupled, the feedback de- 
signed to suppress the leading mode will always affect, 
and can potentially destabilize, the subleading modes. 
As our numerical calculations show, such destabilization 
does indeed occur, unless the spatial profile w{x) of the 
thermal perturbation is chosen carefully to avoid this. 
Had the evolution operator been normal, we could have 
always chosen the feedback in such a way that different 
modes became uncoupled, so only the stability of a few 
independent unstable modes had to be considered. The 
problem of controlling the contact line instability turns 
out to be quite delicate in comparison. 

Assuming that w{x) is chosen such that all subleading 
modes remain stable, we can truncate the system (34) by 
discarding all modes except the leading one. The stability 
of the single- mode truncation 

Go = (^00 + q^Boo + q^Coo)Go (36) 

is a necessary condition for the stability of the full system 
(34). The truncated system is stable when the following 
three conditions are satisfied: 

Aoo = kw'{oo) < 0, 

-Boo = / ho{ho — l)dx — k / N2{ho)wdx<0, 
Jo Jo 

Coo = < 0. (37) 
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-100, 



FIG. 10. The adjoint cigonfunctions and fi{x) for 

q = 0, a — 0.01 and c = 1. Qualitatively similar profiles are 
obtained for other values of the contact angle. A portion of 
the eigenvalues computed on the interval of length = 80 is 
shown. 



The first condition is satisfied for any feedback local- 
ized near the contact Une, because in this case ^oo — 0; 
while the second condition can always be satisfied with 
the proper choice of the gain k. The third condition is 
clearly satisfied as well. Since (34) is valid for all values 
of q for which the governing equation (9) is valid and 
the choice of k in (37) is independent of q, we can im- 
mediately generalize to non-monochromatic disturbances 
by integrating over all q, such that the feedback will be 
given by 



Ae{x, y, t) = -kT{dea)-^w{x + ut)i{y, t), 



(38) 



where ^(r/, t) is the instantaneous deviation of the contact 
line from its mean position. 

How can we choose the heating profile in the flow di- 
rection which will stabilize the unstable mode without 
destabilizing any of the initially stable modes? As Fig. 8 
shows, the leading mode is localized under the capillary 
ridge. Therefore, to suppress the instability one has to 
apply a perturbation which will be similarly localized to 
within the region occupied by the ridge. The inspection 
of matrix elements (35) shows that such a localized per- 
turbation will affect all modes whose adjoint eigenfunc- 
tions are not small in that region. According to Fig. 9 
both /o and /i are of the same order of magnitude there, 
while the other subleading modes are several orders of 
magnitude smaller. As a result, only the stability of the 
leading and the first subleading mode may change in re- 
sponse to feedback. The numerically computed spectra 
of the controlled system support this conclusion. 

The stability of the two leading modes is determined 
by the signs and magnitudes of the matrix elements Anm 
and Bnm with n,m = 0,1. A few general comments 
about these matrix elements can be made based on the 
structure of the eigenfunctions. As Fig. 10 shows, /o and 
/i are both nearly constant under the ridge and have op- 
posite signs, while go{0) and 51(0) have the same sign. 



As a result, a decrease in Bqq is necessarily accompanied 
by a commensurate increase in Bn. Furthermore, since 
/o is constant, we have Aon = for any n, so at q = 
the eigenvalues of the controlled system are P'q{0) = and 
(i'i{0) « (3i{0)+Aii. The matrix element An is generally 
nonzero and changes in response to the strength of the 
applied feedback. An inappropriate choice of the profile 
u'(x) or the gain constant k can make An large enough 
to cause destabilization at long wavelengths. However, 
even when 'w{x) is chosen such that An is negative (but 
small), destabilization of the n = 1 mode at short wave- 
lenghts can occur due to the increase in Bn, if the feed- 
back gain is too large. 

Additional insights can be gained by considering the 
effect of feedback from the physical point of view. The ac- 
tion of feedback in the direction transverse to the flow is 
described by the influence function N2{ho). The first and 
second term of this function describe the motion of the 
liquid under the action of, respectively, surface forces and 
pressure produced by the local gradients in the surface 
tension. For a convex region of the film, where < 0, 
these two effects will compete with each other. For in- 
stance, a local maximum of surface tension will induce 
the flow along the surface toward that location and the 
flow in the interior of the liquid layer away from that 
location. The first effect will dominate for low capillary 
ridges (small contact angles), the second one for high 
capillary ridges (large contact angles). 

As the instability is caused by the increased mobility of 
the capillary ridge, one might envision enforcing control 
by changing the thickness of the film. The local thick- 
ness of the capillary ridge could, in principle, be modified 
by locally heating or cooling it to redistribute the liquid 
in such a way as to decrease the thickness, and hence 
the mobility, where we need to slow down the motion of 
the contact line and increase the thickness and mobility, 
where we need to speed it up to compensate for the de- 
viation. This can be achieved by choosing w{x) which is 
localized under the ridge and does not change sign. For 
instance, one could pick a Gaussian profile representing 
the effect of thermal spreading in the solid substrate 



w{x) = exp 



{x - Xq)^ 
2Aa;2 



(39) 



where xo and Aa; are chosen such that w{x) is centered 
under the capillary ridge and has a comparable width. 
A sample profile is shown in Fig. 11(a) for a special 
choice of parameters. The numerically computed spec- 
trum. Fig. 11(b), demonstrates that the stabilization can 
indeed be achieved by this method for driven films with 
a relatively low capillary ridge. 

This simple approach, however, does not always 
achieve the desired result. In fact, it only succeeds when 
the largest growth rate l3o[qmax) in the uncontrolled sys- 
tem is small compared to |/3i(0)|, i.e., when the feedback 
required to stabilize the n = mode is too weak to desta- 
bilize the n = 1 mode. One can already notice the sign of 
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FIG. 11. Stabilization of a thermally driven film with a low 
capillary ridge, hmax ~ 1-63 {a = 0.1 and c — 1). (a) The 
asymptotic state ho, influence function N2{ho), and thermal 
perturbation profile w. (b) The two leading eigenvalues of the 
original and controlled system for A; = 1. 



approaching trouble by looking at Fig. 11(b): the eigen- 
value of the n = 1 mode at q = starts to creep upward 
due to the increase in the matrix clement An. If this 
approach is used to stabilize a flow with f3o{qmax) com- 
parable to |/3i(0)|, the feedback required to suppress the 
n = mode becomes strong enough to destabilize the 
n = 1 mode at low wavenumbers. The numerically com- 
puted spectra for higher capillary ridges {hmax > 1-7) 
show that the low wavenumbers are destabilized before 
the high wavenumbers are stabilized for any choice of the 
width Ax of the thermal perturbation. 

We are thus forced to look for an alternative solution. 
Changing the overall thickness of the capillary ridge, even 
locally, is a rather ineffective procedure, especially for 
small q, as the liquid has to be redistributed over large 
distances in the transverse direction. One could instead 
apply a local force to the ridge, redistributing the liquid 
between its front and back. For instance, by heating the 
front of the ridge and cooling its back one creates the 
pressure gradient enhancing the flow. Reversing the sign 
of the applied perturbation impedes the flow, directly 
affecting the propagation velocity. The corresponding 
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FIG. 12. Stabilization of a thermally driven film with a 
high capillary ridge, hmax = 2.75 (a = 0.01 and c = 1.95). 
(a) The asymptotic state /lo, influence function N2{ho), and 
thermal perturbation profile w. (b) The two leading eigen- 
values of the original and controlled system for fc = 4. The 
n = 1 mode is strongly suppressed by feedback, the respective 
eigenvalue lies outside the graph. 



thermal perturbation should have a profile 'w{x) which 
changes the sign near the highest point of the ridge. For 
instance, if one chooses the profile to be anti-symmetric 



w{x) 



xo ) exp 



{x - Xq)^ 
2Aa;2 



(40) 



the matrix element An can be made large and negative, 
so the n ~ mode can be stabilized without destabilizing 
the n ~ 1 mode. Fig. 12(a) shows that such a thermal 
perturbation with the position xq and width Ax tuned 
to be roughly the same as those of the capillary ridge, 
should have the largest effect on Bqq. Indeed, the influ- 
ence function N2{ho), dominated for high capillary ridges 
by its second term, also changes sign near the highest 
point of the capillary ridge. The numerically computed 
spectrum. Fig. 12(b) shows that one can again success- 
fully suppress the instability. 

Another alternative is to exploit the narrow concave 
region of the film near the contact line. One can again 
use the Gaussian thermal profile (39) centered at the con- 
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FIG. 13. Stabilization of a thermally driven film with a 
medium height capillary ridge, hmax = 1.93 (q = 0.01 and 
c = 1). (a) The asymptotic state ho, influence function 
N2{ho), and thermal perturbation profile w. (b) The two 
leading eigenvalues of the original and controlled system for 
k = 4. The n = 1 mode is strongly suppressed by feedback, 
the respective eigenvalue lies outside the graph. 



tact line (see Fig. 13(a)). Heating this region and thus 
lowering the surface tension one produces gradients in 
both the pressure and surface tension, which induce the 
secondary flow away from the contact line. Respectively, 
cooling this region draws the liquid towards it, providing 
a direct way to locally control the propagation speed of 
the film. Making the amount of heating or cooling pro- 
portional to the displacement of the contact line again 
allows one to suppress the contact line instability. The 
numerically computed spectrum of the system with and 
without feedback is shown in Fig. 13(b). 

Finally, let us look at the transient amplification of 
disturbances in the presence of feedback. It is unclear a 
priori what effect the control would have on the transient 
dynamics. Numerical calculations show that, depending 
on the choice of the thermal profile w{x), feedback can 
either increase or decrease; the degree of nonnormality in 
the system. This observation is consistent with the the- 
ory. As the matrix elements (35) show, feedback affects 
different modes in a different way, and small changes in 




FIG. 14. The transient amplification with and without con- 
trol for a = 0.01, c = 1.95, fe = 4, and the anti-symmetric 
thermal perturbation profile (40). 



the eigenvalues can have a large effect on the transient 
amplification. For instance, thermal perturbations with 
an anti-symmetric profile, such as (40), can decrease the 
transient amplification at small wavenumbers by almost 
an order of magnitude (see Fig. 14). Generally, direct 
control of the propagation velocity via longitudinal sur- 
face tension gradients leads to a decrease in the transient 
amplification, while indirect control via transverse sur- 
face tension gradients affecting the mobility increases the 
transient amplification for long wavelength disturbances, 
which is also consistent with naive expectations. 

We can thus conclude that an appropriately chosen 
feedback can make the dynamics asymptotically stable 
without increasing the transient amplification of distur- 
bances, such that the contact line instability is sup- 
pressed in the presence of noise characteristic of typical 
experimental conditions. Moreover, the feedback is capa- 
ble of reducing the transient amplification as well, so we 
can expect that feedback control can be effective in sup- 
pressing the contact line instability even when it is caused 
by the nonlinear effects. Finally, once the instability is 
suppressed, selective patterning can be achieved by re- 
moving feedback and/or introducing additional forcing 
at a wavenumber corresponding to a desired pattern. 

The proposed control algorithm has been verified 
experimentally. Although no systematic investigation 
of different thermal perturbation profiles has been at- 
tempted so far, the proof-of-principle experiments have 
shown that by heating the advanced regions of the film 
and cooling the retarded regions one can completely sup- 
press the contact line instability. It has also been shown 
that transversely modulated thermal perturbations ap- 
plied near the contact line can be used to achieve selec- 
tive patterning, producing perfectly periodic patterns of 
rivulets. The details of the experiments will be presented 
in a separate publication (N. Garnier, R. O. Grigoriev, 
and M. F. Schatz, "Optical manipulation of microscale 
fluid flow" , in preparation) . 



11 



VI. DISCUSSION 



To summarize our results, we have determined that 
the cvohition operator governing the dynamics of sponta- 
neous disturbances for thermally driven films, both with 
and without feedback, is significantly nonnormal and can 
transiently amplify those disturbances. The strongest 
transient amplification occurs for the zero wavcnumber 
which does not lead to contact line instability. How- 
ever, even for nonzero wavenumbers transient amplifica- 
tion is unlikely to produce an instability for levels of noise 
characteristic of typical experimental conditions. There- 
fore, linear stability analysis accurately describes both 
the controlled and the uncontrolled dynamics. In con- 
trast, for gravity driven films at small inclination angles 
the transient amplification could be much stronger, with 
the maximum achieved at a nonzero wavenumber^^, pro- 
viding an alternative mechanism for instability. 

We have also shown that the contact line instability in 
thermally driven films can be effectively suppressed by lo- 
cally heating or cooling the liquid behind the contact line. 
Such thermal perturbation can be easily imposed exper- 
imentally with sufficient spatial and temporal resolution 
by radiatively heating the substrate^'^. This approach 
offers significant advantages in controlling the dynamics 
of microflows compared to the one based on chemical 
patterning of the substrate^°'^^. First of all, no prepara- 
tion of the substrate is needed, while the patterns can be 
dynamically reconfigured, offering potential for a signifi- 
cant increase in flexibility. Second, feedback control can 
be used to achieve extremely small feature size, if high 
intensity radiation is used to drive the flow on a thin 
substrate with moderate thermal conductivity^, opening 
up new prospects for microfluidics and microfabrication 
applications. Finally, feedback control provides a unique 
opportunity for studying the dynamics of subdominant 
modes and even unstable states of the system. For in- 
stance, it can be used to experimentally measure the 
growth (or decay) rates of monochromatic disturbances 
with wavenumbers qq ^ Qmax by using a wavenumber de- 
pendent gain to suppress the disturbances which would 
otherwise obscure the dynamics of the mode of interest. 

We also expect that feedback control can be equally ef- 
fective in suppressing nonlinear instabilities such as those 
occiirring in gravity driven spreading at small inclina- 
tion angles. Indeed, we have seen that the profile of 
the thermal perturbations can be tuned to decrease the 
degree of nonnormality. Therefore, by suppressing the 
transient growth feedback can also quench the bootstrap- 
ping mechanism leading to a nonlinear instability. How- 
ever, because small disturbances at the contact line could 
be transiently amplified to produce 0(1) changes in the 
thickness of the capillary ridge^^, it is possible that the 
control algorithm will have to use direct measurements 
of the thickness rather than the much easier to monitor 
position of the contact line. 
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